Non-modal approach to linear theory: 
marginal stability and the dissipation of turbulent fluctuations 
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The non-modal approach for a linearized system differs from a normal mode analysis by following 
the temporal evolution of some perturbed equilibria, and therefore includes transient effects. We 
employ a non-modal approach for studying the stability of a bi-Maxwellian magnetized plasma 
using the Landau fluid model, which we briefly describe. We show that bi-Maxwellian stable 
equilibria can support transient growth of some physical quantities, and we study how these 
transients behave when an equilibrium approaches its marginally stable condition. This is relevant 
to anisotropic plasma, that are often observed in the solar wind with a temperature anisotropy close 
to values that can trigger a kinetic instability. The results obtained with a non-modal approach are 
relevant to a re-examination of the concept of linear marginal stability. Moreover, we discuss the 
topic of the dissipation of turbulent fluctuations, suggesting that the non-modal approach should 
be included in future studies. 
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I. INTRODUCTION 



Linear theory represents a powerful tool for the 
interpretation and understanding of many space plasma 
properties observed by in situ spacecraft. For instance, 
the temperature anisotropy of the solar wind is thought 
to be bounded by values that are consistent with the 
stability thresholds derived within the linear theory 
of the Vlasov-Maxwell set of equations. The broadly 
accepted view is that a macroscopic property of the 
plasma (such as temperature, density, or mean velocity) 
can be constrained by the nonlinear feedback associated 
with a linear instability. This is because the primary 
consequence of any instability is to reduce the amount 
of free energy that drives the instability, relaxing the 
plasma towards a marginally stable condition. This is 
equivalent to saying that the plasma is unlikely to be 
found in an unstable state because it tends to change 
its macroscopic properties in a way that would lead to a 
linearly stable condition. 

In the solar wind it is argued that the expansion of the 
plasma from the Sun in a radially decreasing magnetic 
field should produce much higher values of temperature 
anisotropy (with respect to the background magnetic 
field) than those observed. Also, it has been shown that 
the highest values of anisotropy observed in the solar 
wind are consistent with the thresholds of linear kinetic 
instabilities driven by temperature anisotropies. This is 
true both for protons [H, H3, [2(| and electrons [H, HtJ , 
for a large range of plasma beta, and both for 2j_/Tn > 1 
(whistler and mirror instabilities), and for Tj_/T» < 1 



(firehose instability) . 
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The physical behaviour of an unstable anisotropic 
collisionless plasma subject to the electron firehose 
instability, that relaxes towards marginal stability, has 
been elucidated in Camporeale and Burgess with 
fully non-linear Particle-in-Cell simulations. What 
emerges is that the plasma state is likely to be found 
bouncing around the marginal stability threshold, due to 
the competition of two mechanisms: the reduction of the 
anisotropy and plasma beta caused by the development 
of the firehose instability (above the threshold), and 
the increase of these quantities due to the damping of 
magnetic fluctuations (below the threshold) that result 
in the energization of the particles, predominantly in 
one direction. 

Despite the success of the linear theory to delineate the 
macroscopic properties in which the solar wind plasma 
is more likely to be found (i.e. in stable conditions), and 
the good agreement between linear theory predictions 
and solar wind data, at least two contradictions remain 
unquestioned. 

First, the greater part of the protons and electrons in the 
solar wind, at any distance from the Sun, is observed to 
be isotropic or very lightly anisotropic [lg, [37]], i.e. in a 
condition where no anisotropy instability can be excited 
and therefore the aforementioned argument associated 
with the linear constraining mechanism cannot be 
invoked. In other words, the fact that linear instabilities 
do not allow the plasma to develop anisotropies higher 
than certain values, does not explain why most of the 
plasma is found to be very far from those values, since 
the effect of the expansion should result in a continuous 
increase of the anisotropy. Observational results have 
suggested, as a possible explanation, that the collisional 
age is related to the isotropization of thermal electrons 
[33 . H3|, but the relative importance of collisions and 
instabilities is still unknown. 

Second, a relatively high occurrence of short wavelength 
magnetic fluctuations with small amplitude is persis- 
tently found in the plasma, even in stable conditions. 
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Unless one assumes that a turbulent cascade is able 
to produce such fluctuations at a rate that perfectly 
balances the dissipation, this might be in contradiction 
with the notion that a linear perturbation in a stable 
plasma would damp exponentially in time. 
The aim of this paper is to present a non-modal approach 
to the linear stability problem for a collisionless plasma 
in an uniform magnetic field that, by offering a new 
framework for the understanding of linear marginal 
stability, will reconcile those two apparent observational 
contradictions with a consistent physical interpretation. 
Another important application of the plasma linear 
theory that will be revisited in the light of a non-modal 
approach is the interpretation of the physical mechanism 
that controls the kinetic damping of fluctuations in 
the turbulent dissipation range. This is an area of 
increasingly active research, and it has been argued 
that the linear approximation (i.e. the assumption that 
perturbed quantities have a much smaller amplitude 
than the equilibrium ones) might be used in this scenario 
[To| . There is a general agreement on the fact that the 
nonlinear cascade of turbulent fluctuations should take 
into account the onset of kinetic effects at a certain spa- 
tial scale, around the ion Larmor radius [27j ■ It is indeed 
observed that the power density of magnetic fluctuations 
undergoes an abru pt steep ening for wavenumbers above 
a certain value (l6l. \24. l3l| . and it is thought that kinetic 
effects might be responsible for the steepening. However, 
a complete understanding on exactly why this happens 
and what are the important parameters that determine 
at which scale one should expect the kinetic effects to 
produce a steeper slope of the power spectrum (which 
is observed to follow a power law in wavenumber both 
above and below the steepening) is still missing. 
Altough the study of kinetic instabilities, and their ap- 
plication to the solar wind plasma, has been traditionally 
kept separate from the studies addressing the turbulent 
dissipation range, we look at those two aspects of plasma 
physics as interlocked. In fact, they describe the same 
problem from two distinct viewpoint, because the un- 
derstanding of turbulent dissipation at small scales can 
be properly addressed only by including in the energy 
balance also the injection of turbulent fluctuations due 
to the development of kinetic instabilities. 
Different approaches to t his p roblem include the use of 
gyrokinetic linear theory [ljj an d simulations [2l| . 
Hall-MHD d, S3, ii| and Particle-in-Cell simulations 
[3]. A simple diffusive model for the nonlinear cascade 
in the inertial range (i.e. for spatial scales that can 
be precisely studied with a MHD description) was 
suggested by Zhou and Matthaeus [4 If. This model has 
been successively applied by Li et al. [25| for addressing 
the transition between inertial and dissipation ranges. 
Their conclusion was that any kinetic linear damping 
mechanism would not be able to reproduce the observed 
power spectra of magnetic fluctuations, because it would 
produce a steep cut off instead of a power law, for high 
wavenumbers. 



We will show that, by approaching the linear theory with 
a non-modal formalism, the diffusive model proposed by 
Li et al. (2f| does not necessarily rule out the possibility 
that a completely linear damping mechanism might be 
responsible for the ultimate dissipation of turbulent fluc- 
tuations. Moreover, we argue that the results based on 
a non-modal approach clarify the relationship between 
different eigenmodes, and clearly establish the fact that 
the problem of turbulent dissipation is very unlikely to 
be understood by studying the damping properties of 
one single normal mode, as previous works attempted to 
do. We will argue that, at kinetic scales, the evolution 
of a small perturbation depends very much on the linear 
interactions between many different modes, including 
heavily damped ones. Also, from a purely computational 
point of view, we will show that our results question the 
possibility of exciting one single eigenmode. 



A. Modal and non-modal approaches 

Plasma stability theory is traditionally treated as 
a normal mode ('modal') analysis. This means that 
the focus is on the eigenmodes of a slightly perturbed 
equilibrium, which are assumed to grow or damp 
exponentially in time. The most rapidly growing (or 
the least damped) eigensolutions are the object of 
greatest interest in the normal mode analysis, which 
therefore investigates the stability problem purely in its 
time-asymptotic solutions, regardless of the details of 
the perturbation imposed on the initial equilibrium. 
The non-modal approach differs from the normal mode 
analysis in treating the linear stability as an initial value 
problem. This allows the study of transient phenomena 
and to follow the evolution of a perturbed system in 
time. This evolution depends, of course, on how the 
initial equilibrium is perturbed. 

The non-modal approach to the plasma stability problem 
can in some cases be crucial to the understanding of the 
evolution of a plasma in the linear regime, and it should 
therefore routinely be used. Especially when studying 
a system in stable conditions the normal-mode analysis 
misses the phenomenon of transient growth in time of 
some physical quantities. This is a well-known effect 
which has been long studied in hydrodynamic flows 
[33l ]. It is related to the spectral properties of the linear 
operator that describes the evolution of the system in 
time. In particular, when an operator is non-normal, 
i.e. it does not commute with its adjoint, the norm of 
an initial perturbation applied to the equilibrium can 
grow in time by large factors, before decaying, even if 
all the eigenmodes of the operator are damping (i.e. the 
equilibrium is stable). 

Despite the fact that transient growth has been known 
of for a long time , it has not been emphasised in the 
theory of plasma stability. We have conducted the first 
study of transient growth for a stable kinetic plasma 
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in a previous paper showing that this phenomenon 
appears to be related to the kinetic regime of a plasma: 
it is indeed more accentuated for higher plasma beta and 
shorter spatial scales. In order to include kinetic effects 
we used a Landau fluid model, which incorporates linear 
Landau damping, and finite Larmor radius corrections. 



B. Aims of the paper 

The purpose of this paper is threefold. First, we ex- 
tend the work of Camporeale et al. [5( to an anisotropic 
plasma, for conditions typical of the solar wind. By doing 
so we will highlight the inadequacy of the normal mode 
analysis for a stable kinetic plasma. 
Second, we will clarify the relation between transient 
growth and marginal stability, and we will propose a 
physical picture where the two observational contradic- 
tions mentioned above will be reconciled with the linear 
theory. 

Third, we will discuss the importance and appropriate- 
ness of a non-modal linear theory for the understanding 
of the dissipation of turbulent fluctuations. It should be 
the goal of a future complete theory of solar wind tur- 
bulence to elaborate a unified theoretical framework that 
links the physics of kinetic instabilities, and the damping 
of turbulent fluctuations in a consistent manner, and we 
believe that such a theory will benefit from embracing a 
non-modal approach to the stability problem, such as the 
one presented in this paper. 
The paper is organised as follows. 

In Section 2 we describe briefly the Landau fluid model, 
that has been used to obtain the results presented in the 
paper. The model has been described and commented 
in length elsewhere [28|, |29|, |38| , hence we will just briefly 
present the main features of the model, and we report the 
set of equations in the Appendix. In Section 3 we will in- 
troduce the methodology of the non- modal approach, and 
we will define some important quantities for our analy- 
sis. Section 4 describes the results of our study, with 
emphasis on the relationship between transient growth 
and marginal stability, and on the non-modal approach 
for the study of turbulent dissipative fluctuations. A final 
discussion, with suggestions for future work, is reported 
in Section 5. 



II. THE FLR-LANDAU FLUID MODEL 

The idea of incorporating Landau damping in a set of 
fluid equations was introduced by Hammett_and Perkins 
(l7| and later developed in Snyder et al. [36| . where the 
fluid hierarchy obtained from the drift kinetic equation 
is closed at the level of the third or fourth order moment. 
These models are limited to scales large compared with 
the ion gyroradius. Other models, called gyrofluid mod- 
els, consider the fluid hierarchy obtained from the gy- 



rokinetic equation, providing a set of equations for fluid 
moments suitable for the description of sub-Larmor ra- 
dius scales [l|, H, 0]. The equations of gyrofluid models, 
however, are not written in the physical coordinates but 
in the gyrocenter variables, making their interpretation 
more difficult. A simpler formulation retaining hydrody- 
namic nonlinearities together with a linear approxima- 
tion of FLR contributions was recently developed by de- 
riving equations for the hydrodynamic moments directly 
from the Vlasov-Maxwell system [H [H, HI- This is 
the model used in this paper. The hierarchy of fluid 
equations is closed at the level of the fourth order mo- 
ment. In its linearized version, the plasma dispersion re- 
lation is approximated by a suitable Pade approximant, 
and this allows to cast the linear set of equations as 
a standard eigenvalue problem. In addition to the hi- 
erarchy closure, this approach involves the modeling of 
FLR effects by expressing the non-gyrotropic part of ten- 
sors such as pressures, heat fluxes, or fourth order mo- 
ments in terms of lower-rank moments, in a way consis- 
tent with the linear kinetic theory in the low-frequency 
limit e ~ oj/VLi <C 1, for both quasi-transverse fluctu- 
ations (k\\/k± ~ e) with no condition on k±r^ (as in 
gyrokinetic and gyrofluid approaches), but also for hy- 
drodynamic scales with hi ~ k± <C 1/V/,. Here tti de- 
notes the ion cyclotron frequency and the ion Larmor 
radius. At large scales, the model, which then reduces 
to usual anisotropic MHD, also captures the fast waves, 
in contrast with gyrofluids. The frequency and damping 
rates of low-frequency waves are accurately described in a 
range of scales that extends to small (sub-Larmor radius) 
scales when the propagation direction is almost perpen- 
dicular to the ambient magnetic field (according to the 
gyrokinetic scaling). The complete model is quite in- 
volved and is thoroughly described in Passot and Sulem 
[29| . The full set of equations is reported in Appendix, 
for completeness. The total number of variables is 16, 
hence the linear problem reduces to the formulation of a 
16 x 16 complex matrix. 



III. NON-MODAL APPROACH 

In this section we introduce the mathematical tools 
and the methodology that will be used in the rest of the 
paper. For a more complete description of the non-modal 
stability theory (in the context of hydrodynamics) we re- 
fer to the review by Schmid [33| : several issues related 
to the non-normality of a linear operator (that will be a 
central point of our discussion) have been described in 
gre at depth in the monograph by Trefethen and Embree 

The Landau fluid (LF) model can be linearized as usual, 
by writing each physical quantity as a sum of an equi- 
librium and a perturbed contribution (subscript and 1 
respectively): 4> = 4>q + e<pi, assuming that e < 1, and 
neglecting terms of order higher than one in e. Once 
the first-order variables are Fourier decomposed (drop- 
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ping the subscript) <p(r, t) = 0exp[i(k • r)], the linear LF 
model can be cast as a set of ordinary differential equa- 
tions for the complex amplitudes 4>: 



(1) 



A is an autonomous (i.e. not a function of time) op- 
erator, that takes the form of a 16x16 sparse complex 
matrix, whose entries depend on the properties of the 
plasma (protons and electrons temperature anisotropy 
and plasma beta), and on the magnitude and angle of 
propagation of the wavevector k. 

The fact that the LF model describes the plasma through 
a set of only 16 variables makes the analysis of the ma- 
trix A computationally affordable without any particular 
method used for large matrix manipulations. Accord- 
ingly, all the results presented in this paper have been 
produced using in-built routines of MATLAB. 
The solution of Eq. |T]) is given by 



0(t) = e A V(0) 



(2) 



where </>(0) is the state vector of the initial perturbation 
and the exponential of the matrix (which is defined as 



Ai 



At + \{Atf 



.) completely determines 



the evolution in time of the initial state (it is sometime 
called the propagator). 

A distinctive feature of the matrix A that is crucial to 
our argument is its normality. If A commutes with its 
adjoint AA* = A*A, than A is said to be normal, and 
its eigenvectors form a complete orthogonal set. We 
note that the definition of normality depends on which 
norm one refers to. Any non-normal linear operator can 
be made normal by choosing a different definition of the 
norm, from which the definition of adjoint follows Q. 

We use here the 2-norm ||u|| = \Jj2i=i l u «l 2 - Altough 
this norm does not provide an information on single 
variables, it indicates whether the perturbation complies 
or not to the requirement of being much smaller than 
the equilibrium quantities. What is important is that a 
large norm of a perturbation implies a deviation from 
the assumption of linearity, and therefore might result 
in the triggering of non-linear effects. 



A. Spectral abscissa and numerical abscissa 

In order to study how the plasma responds to a small 
perturbation we introduce the quantity G(t), which mea- 
sures the amplification (or reduction) in time of an initial 
perturbation: 



G(t) 



\m)\\ _ lle A V(0) 
II0(O)|| 11^(0)11 



(3) 



This quantity clearly depends on the details of the initial 
perturbation </>(0), but it is always bounded by above 



from the quantity ||e A *||, which defines the supremum 
of G(t) for any possible 0(0). The behaviour of ||e A *|| 
depends on the spectral properties of A. If one indicates 
with er(A) the spectrum of A, i.e. the set of z G C such 
that (zI — A) is singular, where I is the identity operator, 
then the quantity 



a(A) = maxSR[cr(A)] 



(4) 



is referred to as the spectral abscissa. To determine this 
quantity is the main objective of the normal-mode analy- 
sis, because it gives the growth rate of the most unstable 
mode (for an instability) , or the damping rate of the least 
damped mode (for a stable plasma). Hence, the time- 
asymptotic evolution of any initial perturbation is given 
by e"^- 1 *. Moreover, for a normal operator, the spectral 
abscissa bounds the quantity G(t), for any time t > 0: 

sup G(t) = ||e At || = e Q(A)t (for t > iff A is normal). 

In general however, for a non-normal operator, Eq.^ 
holds only in the limit t — > oo, and there is no exact 
formula for the quantity ||e ||, that can only be approx- 
imately estimated for t > (differe nt ap proximations can 
be found in Trefethen and Embree [391]). 
But a rigorous formula exists for the growth of ||e At || at 
t = 0. This is referred to as the numerical abscissa: 



1(A) = 



super 



{=0 



(6) 



Note that since (A + A*) is Hermitian, 77(A) is real. 
The numerical abscissa provides the information about 
the highest possible growth rate of any initial perturba- 
tion, at time t — 0. For a normal operator it follows from 
Eq.© that 77(A) = a(A), consistently with the fact that, 
if A is normal, the quantity G(t) is bounded by e™^* 
for any time. 

What is surprising is that, for a non-normal operator, 
even if all the normal modes are damping (i.e. a(A) < 0) 
the numerical abscissa can be positive, hence allowing the 
quantity G(t) to grow for some particular initial pertur- 
bations. It has been shown in Camporeale et al. [5| that 
G(t) can indeed reach values of about 10 3 — 10 4 over 
short periods of time, for a Maxwellian plasma described 
by the Landau fluid model. 

The fact that a perturbation could grow while the eigen- 
modes of the operator decay in time appears counter- 
intuitive. It is purely due to the superposition of non- 
orthogonal eigenvectors, and we refer to Figure 2 in 
Schmid (33j for a graphical paradigmatic explanation of 
this effect. 



B. Pseudospectra 

A key aspect of non-normal operators is that the spec- 
trum may be highly sensitive to small perturbations. In 
order to quantify the effect of small perturbations on a 
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linear operator, the concept of pseudospectra has been 
introduced. There are at least four definitions of pseu- 
dospectra, that have been shown to be mathematically 
equivalent [3{| . We report here the two that have a most 
immediate link with a physical interpretation. If one em- 
ploys the convention that the spectrum <r(A) (the set 
of eigenvalues) is formed by complex values z for which 
|| (zl — A) _1 || = oo, the e-pseudospectrum cr £ (A) is de- 
fined as the set of z 6 C such that 

\\(zl - A)" 1 1| > - (Definition 1) 

£ 

for any e > 0. 

This is equivalent to saying that z is an eigenvalue of the 
perturbed operator (A + E) for some operator E with 
||E|| <e, i.e.: 

z e cr e (A) ^=^> z e er(A+E) with |[E|| < e (Definition 2) 

Let us now comment on these definitions and their mean- 
ing in a more physical sense. Our linear operator stud- 
ies the evolution of a small perturbation applied to an 
equilibrium, namely a bi-Maxwellian plasma in an uni- 
form magnetic field. However, one could argue that 
small discontinuities or inhomogeneities of the equilib- 
rium quantities might be modelled as small disturbances 
to the linear operator. The question that pseudospec- 
tra quantitatively answer is: how does the set of eigen- 
values change under the effect of small perturbations of 
the operator ? From definitions 1 and 2, one see that 
the e-pseudospectrum a £ tends precisely to the standard 
spectrum a, when e — > 0. In the complex plane the e- 
pseudospectra are the open subset that contains all the 
eigenvalues of the perturbed operator (A + E), for any 
possible perturbation E, such that ||E|| < e. Therefore 
they give a measure of the distortion of the spectrum 
due to the perturbation applied to the operator. It is 
straightforward to prove that the e-pseudospectra are a 
nested set. That is, <r El (A) C er £2 (A) for e% < e-i 
An example of pseudospectra of our LF operator is 
given in Figure [TJ for an isotropic plasma with = 10 
and k = 0.5. Each contour corresponds to values of 
e = 10~ 6 - 8 , 10~ 6 - 6 , • • • , 10" 5 - 6 . Only a part of the com- 
plex plane is shown, with 9 of the 16 eigenvalues visible. 
The important point here is to understand that pertur- 
bations of a certain size make the e-pseudospectrum so 
large, that it would contain many eigenvalues of the orig- 
inal unperturbed operator. In a sense, this means that 
such perturbations distort the problem in such a way that 
the information about the modes of the unperturbed op- 
erator becomes useless. In this respect the different be- 
haviour between normal and non-normal operator is most 
evident. The e-pseudospectrum of a normal operator is 
defined as the union of open balls about the points of the 
spectrum: 



of the pseudospectra can be computed straightforwardly 
once the eigenvalues are known. We show in Figure [TJ 
with dotted lines, how the contours of a £ would be for 
the same LF operator, if it were normal, for e = 0.00075. 
This contour is qualitatively similar to the (true) contour 
for e = 10~ 5 ' 8 , embracing all the 9 eigenvalues. This is 
the crucial point about the different behaviour between 
normal and non-normal operators. A perturbation of the 
order 10 -5 8 is sufficient to achieve a distortion of the 
spectrum that, if the operator would have been normal, 
would have required a perturbation of the order 7.5-10~ 4 , 
i.e. about 470 times higher. 

The displacement of the spectrum of an operator sub- 
ject to small perturbations implies two facts. On one 
hand some, if not all, of the modes become somehow 
coupled, i.e, their relative distance in the complex plane 
can be modified, and their properties (like phase speed 
and damping rate) changed. On the other hand, from 
a computational point of view, the non-normality of the 
operator almost completely rules out the possibility of 
exactly exciting one single eigenvector. Small errors, due 
for instance to digits truncation or approximation, can 
result in the excitation of a 'pseudomode', that could lie 
in the complex space, very far from the mode that was 
intended to be excited, and that is in reality a super- 
position of different non-orthogonal modes. This might 
lead to an initial transient behaviour, which is not un- 
usual in numerical simulations, even though it is seldom 
commented. 



IV. RESULTS 

In this section we apply the linear LF model to a 
stable proton-electron (p, e) bi-Maxwellian plasma in an 
uniform magnetic field. We are interested in how the 
distance from marginal stability affects the evolution of 
linear perturbations in the plasma, and wish to address 
transient behavior that is missed by the standard modal 
analysis. For all the results shown in this section the 
protons are considered isotropic with T± p = Tj| p = Tj| e . 
All quantities are therefore referred to electrons, and the 
subscript e is dropped. 

We focus on quasi-perpendicular waves, since the LF 
model has a domain of validity that extends to small 
scales only for oblique wavevectors. The angle of prop- 
agation for the wavevector k with respect to the back- 
ground magnetic field is 9 = tan _1 (1000). 
It is known that the marginal stability thresholds (i.e. the 
curves for which the growth rate is exactly null) obey a 
law of the form: 5r=- = 1 + -J7, where S and a are con- 

T\\ p 

stants. We have derived the stability thresholds for the 
LF model, and we have computed the fitting parameters 
using a Levenberg-Marquardat method (30j . The result 
is: 



where dist indicates the distance of a point to a set in 
the complex plane. For a normal operator the curves 
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and 

^ = 1 + ^ (far II, < 31) (8) 
-Ml P 

The values for Tj| > T± are not dissimilar from those 
obtained from the Vlasov-Maxwell equations for the 
electron firehose instability (S = —1.29, a = 0.98), 
which is known to yield highest growth rate for quasi- 
perpendicular propagation Q. The parameters for 
Tii < Tj_ are instead not comparable with those obtained 
for the mirror and whistler instabilities [121 ]. because in 
this case the whistler instability, which is parallel propa- 
gating, sets the instability threshold [ll[ , and the Landau 
fluid model becomes invalid for strictly parallel wavevec- 
tors. We therefore will focus mainly on the case Ty > T±. 



A. Transient growth and marginal stability 

It has been shown in Camporeale et al. Q that an 
isotropic Maxwellian plasma is able to sustain large tran- 
sient growth of an initial perturbation, that will eventu- 
ally decay in later times. 

As we have seen in Section 3, the spectral abscissa a 
(which is negative, for a stable plasma) provides the in- 
formation about the late-time damping rate of the per- 
turbation. That is, any initial fluctuation will decay as 
e Qt , for large times. An interesting point is to actually 
analyse the time at which the system starts to behave as 
predicted by the normal mode-analysis. 
We show in Figure [5] the amplitude of the y component 
of the magnetic field, normalized to its initial value, for 
two particular initial conditions. The wavevector k is 
chosen equal to 1, and T_i_/Tj| — 1. The parameter [3 is 
equal to 1 in the top panel, and (3 = 5 in the bottom 
panel. What emerges is that, for both cases, the behav- 
ior is highly oscillatory, and no damping is evident before 
TVli — 10 5 . This implies that, for the initial perturba- 
tion to decay, the plasma, once perturbed, should evolve 
without encountering any further perturbation for a very 
large time, which is quite unrealistic. The damping pre- 
dicted by modal analysis (i.e. e at ) is shown in dashed 
line. 

Before proceeding, let us clarify how the initial pertur- 
bations of Figure [5] were chosen. Let us recall that the 
singular value decomposition (SVD) of an operator A is 
given by 

A = U£V*, (9) 

where U and V are unitary matrices, and S is a diagonal 
matrix, that contains the singular values of A. 
If one wants to find the initial perturbation 0(0), that at 
a specific time O undergoes an amplification G(Q) equal 
to ||e Ae || (i.e. the highest possible amplification at time 
0), it is sufficient to calculate the SVD of e Ae , and to 
identify the column vector of V associated with the high- 
est singular value. The initial perturbations chosen to 



produce Figures [2] and [3] are the ones that approximately 
attain the highest amplification at a certain time. One 
could argue that these initial conditions are very special, 
and cannot represent the totality of all the possible initial 
states. This is certainly true. However, it is not feasi- 
ble to apply the non-modal approach to a very large set 
of different conditions (that will never be large enough 
to represent 'all' cases). Beside, the methodology that 
we use, i.e., to focus on only certain particular cases that 
represent the 'worst case scenario', is common also to the 
modal approach, whose results are supposed to be valid 
for any initial conditions (albeit restricted to asymptotic 
times), but still represent only the scenario in which the 
fastest growing (or the least damped) mode is actually 
excited. Later we describe an attempt to approach the 
problem in a more statistical manner. 
Figure [2] shows the results of two cases for an isotropic 
plasma, with different [3, but what is perhaps more in- 
teresting is to look at what happens when the plasma is 
closer to marginal stability, given by Eqs. ([7]) and (jHJ). 
For (3 = 5, Eq. (O predicts that the plasma is marginally 
stable when T±/Tu ~ 0.6. We show in Figure [3] three 
cases of transient growth with (3 = 5 and T±/Tn = 0.65, 
for k = 1,5, 10. The spectral abscissa (damping rate) is 
respectively a = -1.02 • 10~ 5 , -7.2 • 10~ 7 , -3.9 • 10~ 7 . 
Since the damping rate is so low, what is expected is that 
any perturbation applied to the plasma would remain un- 
changed in amplitude for long times. This is indeed what 
is observed. What is rather interesting, however, is that 
the initial perturbations shown in Figure [3] undergo a 
transient growth at early times, and therefore remain am- 
plified for long times. In other words, a transient growth 
effect is able to amplify a perturbation even at marginal 
stability condition, and the fact that the damping rate 
is very small allows to the amplified perturbation to sur- 
vive for long times. In the examples of Figurc[3]thc initial 
value of By gets amplified by a factor between 10 and 100 
and remains so at least for 10 5 ion gyroperiods. 
We now try to quantify the importance of transient 
growth effects in a statistical sense for the following range 
of parameters: (3 £ [1.5, 10], Tx/Ty € [0.2,1.2]. 
We have divided the (/3,T X /I||) space in a 100 x 100 
grid. For each point in the stable region we have gener- 
ated 10000 random initial perturbations. In Figure 0] we 
show the mean value of Git) at times t = 1, 10, 100, 1000, 
in log 10 scale. We recall that G(t) measures the ampli- 
fications of the 2-norm on the vector, hence considering 
all the 16 variables equally. One can notice that at time 
t = 1, Git) increases r for higher /?, almost indepen- 
dently of the value of the anisotropy. Transient growth 
are therefore not a sporadic event, even for an isotropic 
plasma. For later times G(t) has still large values, and 
the highest peaks tend to be localized at the edge of 
marginal stability, for increasing time. In figure O we 
show the average (left panels) and the maximum (right 
panels) value of SB/Bq calculated across the 10000 dif- 
ferent initial perturbations. The highest magnetic fluctu- 
ations are evident in proximity of the Tii > Tj_ instability 
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threshold, for any time. Once again the transient growth 
of magnetic fluctuations appears not to be a sporadic 
and unusual event that happens only for carefully chosen 
initial perturbations, but a feature that is persistent in 
time. 

In the light of these results we suggest the following sce- 
nario that might explain both the presence of fluctuations 
at high wavenumbers in stable conditions, and the fact 
that most of the solar wind is observed close to temper- 
ature isotropy, despite the expansion. 
The increase of anisotropy is expected, in fluid theories, 
due to the conservation of one or more adiabatic invari- 
ants. For instance, the CGL approximation Q predicts 
that the ratio T\\/T± increases as r 2 (where r is the dis- 
tance from the Sun), if the plasma is expanding radially 
in a magnetic field that varies as r -2 . Given that a tem- 
perature anisotropy is a source of free energy for insta- 
bilities, it has been argued that the rise of anisotropy 
should be possible only until the free energy becomes 
so large that an instability might be triggered, with the 
consequence of not allowing a further enhancement of 
anisotropy, and therefore to bound also the increase of 
free energy. 

However, we have shown that, in stable conditions, small 
perturbations can give rise to transient growth of mag- 
netic fluctuations and of high order moments of the par- 
ticle distribution function. We argue that the presence 
of these fluctuations, even far from the stability thresh- 
old, can influence the ability of the plasma to further 
increase the free energy via the anisotropization of the 
particle distribution function. In fact, the thresholds for 
marginal stability are derived from the linear theory as 
the conditions for which the time asymptotic growth rate 
a becomes null, assuming for the plasma to be in equi- 
librium in an uniform magnetic field. No consideration 
is made for the fact that the growth rate is effectively a 
function of time, and that the magnetic field might not 
be uniform, due to transient fluctuations. What might 
happen is therefore that the increase of anisotropy might 
be bounded even before reaching what is considered the 
threshold for anisotropy instabilities. In this scenario, a 
parcel of plasma could experience a 'local' marginal sta- 
bility condition due to a temporary enhanced magnetic 
fluctuation. 

This scenario reconciles the observational contradictions 
mentioned above with a consistent interpretation of the 
linear theory. Of course, trying to interpret space obser- 
vations purely on the base of linear theory might appear a 
naive approach and a complete and coherent understand- 
ing of the process of the solar wind expansion is feasible 
only via a non- linear treatment, and through computer 
simulations. Howewer, our results show that a non-modal 
linear theory could already constitute a strong theoretical 
ground for data interpretation. 



B. The dissipation of turbulent fluctuation 

In this section wc highlight the particular behavior of 
the numerical abscissa as a function of wavenumber, and 
point out certain unique features which indicate that non- 
modal effects may play an important part in the cascade 
of turbulent fluctuations at short wavelengths. 
An approach that is very often used in modelling the dis- 
sipation at short wavelength of turbulent fluctuations is 
to assume that such small fluctuations can be regarded 
as an ensemble of linear waves, each of them damping 
according to their linear damping rate. Howes et al. [201 ] 
have devised a gyrokinetic formalism to study turbulence 
in a magnetized plasma, in the framework of the Goldre- 
ich and Sridhar critical balance assumption [3a |. They 
have shown, by performing nonlinear simulations, that, 
in the case of small damping rates, the linear damping 
does not underestimate the rate at which electromagnetic 
energy is dissipated [2l[. However, the gyrokinetic cas- 
cade model produces an exponential roll-off of the power 
spectrum for high wavenumbers, instead of the observed 
power law [24j |. The fact that a linear damping mech- 
anism cannot reproduce a power-law in the spectrum, 
for high wavenumbers, was previously pointed out by 
Li et al. 25j, using the full Vlasov-Maxwell linear the- 
ory. Our contribution here is to revisit the conclusions 
reached by Li et al. [25] , by showing that a non-modal ap- 
proach might be more appropriate to address this prob- 
lem, which is still largelyunsolved. 

The equation used by [25| to model the spectral energy 
W(k) in wavenumber space is of the kind: 



dW(k) 
dt 



D(k)W{k)+j(k)W{k) + S(k), (10) 



where D(k) is a diffusion operator (that can be modeled 
in such a way to recover a Kolmogorov cascade for large 
wavelengths), S(k) = S(k )S(k — k ) is a source term that 
operates at the (small) injection wavenumber fcoj an d 
7(fe) is the damping rate given by the Vlasov-Maxwell 
linear theory. This model is clearly oversimplified in at 
least three ways. First, it assumes that the nonlinear 
cascade is a diffusive process in wavenumber space. Sec- 
ond, it assumes that the cascade is isotropic in k (while 
it is known that anisotropy plays a very important role) . 
Third, it describes the damping process entirely via an 
ensemble of independent linear waves, all propagating in 
the same direction, and each damping according to its 
respective linear damping rate. 

Despite the weaknesses of the model, it still represents 
a good starting point for our purpose of highlighting the 
importance of transient, non-asymptotic dynamics. The 
main conclusion of Li et al. [25j was that 'a power-law 
spectrum in W(k), in the presence of damping, requires 
j(k) to be a specific power-law'. However, at the light 
of the considerations hitherto drawn, we point out that 
treating the cascade process as a time-asymptotic prob- 
lem rather than an initial value problem, hugely distorts 
the physics. In fact, each wavenumber is continuously 
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subject to injection of energy via the non-linear cascade, 
and to treat the evolution of a fluctuation, as if it was 
injected in a remote past, and let free to damp undis- 
turbed is, in our view, an oversimplification of the whole 
process. 

Although it has been shown that some fluid models are 
able to recover the steepening of the power spectrum 
purely by a nonlinear mechanism (see e.g. Galtier and 
Buchlin Q), we note that if one assumes for the non- 
linear terms in the cascade to act smoothly as a function 
of k, then it is very likely that at least one of the param- 
eters that control the damping at high frequencies must 
present some sort of abrupt change in their dependence 
on k. 

This issue necessitates a deeper and separate study, and 
we do not intend to suggest here a definitive answer. 
However, it is interesting to show the behaviour of the 
numerical abscissa, as a function of k, for different values 
of f3, and T±/Tn (Figure^. We recall that the numerical 
abscissa measures the derivative of the highest amplifi- 
cation that an initial perturbation can undergo, at time 
t = 0. 

Two features are extremely surprising, in Figure[Hl First, 
the numerical abscissa 77 follows a power-law in wavenum- 
bers (i.e. a straight line in the log- log plot). Sec- 
ond, it presents an abrupt steepening of its slope, when 
T±/T\\ 7^ f . As far as we know, this is the only quantity 
derived from the linear theory that presents this two fea- 
tures. Therefore, it makes sense to suggest that the nu- 
merical abscissa should be included in the set of parame- 
ters that control the cascade of turbulent fluctuations at 
short wavelengths. In some sense, this strengthen our ar- 
gument that a non-modal approach to the problem might 
be decisive, and we leave further considerations for the 
future. 



V. DISCUSSION 

Plasma linear theory is traditionally carried out as a 
normal mode analysis. The modal approach, although 
useful to define the conditions for which an equilibrium 
might be unstable, is, in general, unable to describe the 
time evolution of some initial perturbations. In partic- 
ular, it is known that some linear operators can present 
transient growth of physical quantities, even in stable 
conditions. We have addressed the stability of a kinetic 
plasma modelled through a Landau fluid model, via a 
non-modal approach. The following results have been 
established. 

• Some small initial perturbations can undergo large 
amplifications, even if applied to a stable equilib- 
rium (Figure |2J). If the plasma is close to marginal 
stability, this might result in a persistent presence 
of such fluctuations (Figure [3]) . This is because the 
(time-asymptotic) damping rate is so small, that 
the decay of an initial perturbation is appreciable 



only after a time of the order of ~ 10 6 ion gyrope- 
riods. On the other hand, it is quite unrealistic for 
a plasma to be undisturbed for such large times. 

• Transient growth of magnetic fluctuations, or of 
high order moments of the distribution function are 
quite possible also for nearly isotropic plasma, i.e 
far from the thresholds of linear anisotropy insta- 
bilities (Figures 0] and [5} . 

• The numerical abscissa, which is the highest growth 
rate, at time t = 0, for any possible perturbation 
turns out to follow a power-law in k, and to present 
an abrupt change in slope for high wavenumbers. 

We believe that these results could help to interpret 
some space observations in a more consistent way, than 
the traditional modal theory is able to do. We suggest 
that marginal fluctuations observed in stable conditions 
might be the result of transient growth. We also suggest 
that the marginal stability concept should be revisited 
within a non- modal context, and that the prediction of 
the evolution of small fluctuations in a turbulent plasma 
purely by their time-asymptotic damping rate, is an 
oversimplification of the physics. 

Most of the solar wind protons and electrons are observed 
to be nearly isotropic. This is in contradiction with fluid 
theories that predict a continuous anisotropization of 
the distribution function. We suggest that the increase 
in anisotropy might be inhibited by the presence of 
ongoing transient effects, that could limit the capability 
of the plasma to increase the free energy. 
As for the dissipation of turbulent fluctuations, it is now 
acknowledged that the understanding of this complex 
phenomenon, lays in the realm of kinetic plasma physics. 
However, once again the modal approach does not seem 
to be very useful here, and might be misleading. On the 
other hand, the power-law dependence of the numerical 
abscissa, and more importantly the presence of a steep- 
ening for high wavenumbers, is very interesting. We have 
pointed out that, since this is the only quantity (derived 
within linear theory) that presents such features, it 
should be taken into account. 



A. Future directions 

The primary intent of this paper has been to show that 
a non-modal approach could be helpful in reconciling 
some observational evidence with plasma linear theory, 
and to suggest a more frequent use of such an approach. 
Historically, the issues that can arise for a non-normal 
linear operator have been studied in depth in hydrody- 
namics, but have not captured the same attention in the 
space plasma community. 

This of course does not exclude the use of non-linear the- 
ory and computer simulations, that will ultimately be 
the only way to fully understand the processes in action. 
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However, we have shown that the interpretation of simu- 
lations should be made with the awareness of the effects 
produced by non-normal linear operators. 
Altough the Landau fluid model used in this paper is su- 
perior to MHD, being able to capture some kinetic effects, 
it will be desirable in the future to develop a non-modal 
treatment of the full Vlasov-Maxwell set of equations. 
Also, as we have shown, the non-modal approach should 
be framed in a more statistical theoretical framework. 
In conclusion, we believe that the understanding and the 
ability to address transient behaviour will be a crucial 
ingredient in future modelling of space plasma. 



fluxes and gi and the mixed fourth order cumulants 
rlj_. They read 
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APPENDIX A: FLR-LANDAU FLUID SET OF 
EQUATIONS 

In the context of this paper, it is useful to specify the 
form of the model in one space dimension, assuming that 
all the fields only depend on a coordinate £, along a direc- 
tion of the (x, z)-plane making an angle a with the z-axis 
defined by the uniform ambient magnetic field (of magni- 
tude B ). The total plasma density field is normalized by 
po, the magnetic field by B , the velocities by the Alfven 
velocity va — B / '(47rp ) 1 ^ 2 , the pressures by the parallel 
ion pressure po — p l S ^ , the heat fluxes by pqva and the 

fourth rank moments by Pqv a . The unit of length is the 
ion inertial length uaM, and time is normalized to ion 
gyroperiods. The parameter j3 = 8irp /B 2 measures the 
ratio of the (parallel) thermal to the magnetic pressure. 
Velocities without superscripts refer to the ion velocity. 
The electron velocity u e is given by 



u rr = u x + cos a 



d^b y 



U„ = Uy — 



u v = u z — sma- 



(Al) 
(A2) 
(A3) 



d t p + df:{pu) = 
d t (p Uj ) + d s (F} + Ff) = 
d t b p = d^Ey 
d t b y = -d s E p 



dtp\ = -d^{up\ + 



cos a 



|6| 



-q\\ 



+ sinaS| r ) 



+2q r ± V -b - 2p$ • W-6 

cos a „ . „ x 



\b\ 



-q r ± + sin aSjr r ) 



%q\ 



-q r ± W -b-p r ± V -u + p^b-Vu 7 ' -b 

cos a 3/3 r cosa P\ 

-dduq {l + -Trrr u ) - — %(-) 



\b\ 'mi' 2 
-3gfi6- Vu r ■6 + 3fl L V-& 



\b\ 



(A4) 
(A5) 
(A6) 
(A7) 



(A8) 
(A9) 
(A10) 



cos a. (3 r cosa p r ± 
d t q ± = -d,(uq ± + T r ||± ) - -p^—d t {-) 

-q r ± V • u - (^(p| - pi) + r\ {± r~i ± )V • b 



(All) 



(d t + ud s )f\ ± =F- l (- T\T^ l±2 ikq\ - V^Uil*^ 
+f l (f 1 _(T l 1 _ - T^|)cosasina7e|| ± 3 k 2 (d + 1)6^12) 

(dt + «0 £ )rf_ L = F-^-Tp-faikfr - ^{Hl ±1 \k^ ]± 
-TyTl(Tl -Tjj)cosasinaTC| ± 3 k 2 b y ^. (A13) 



where we define b p = cos a b x — sin a b z . We also define 
u = sin a u x + cos a u z which, together with V ■ u = 
d^u, take the same form when using the electron velocity. 
When integrated, the divcrgcnceless condition V • b = 
rewrites cos a b z + sin a b x — cos a, which allows to write 
V-6 = V- (6/|6|) = -cosa<%|6|/|6| 2 . 

The model involves dynamical equations for the ion 
density p and velocity u, the magnetic field components 
bp and b y , together with, for each species r, the gyrotropic 



parallel and perpendicular pressures pjj and p r ± , the heat The electric field and flux terms entering the above equa- 
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tions are given by The operator T 1 is the inverse Fourier transform. The 

second terms on the right-hand-side of Eqs. (|A12p and 

E p = cos a E x - sin a E z (A14) (|A13|1 involve a Hilbert transform with respect to the z- 

_ _/ x ^\ _ }_g p2 (A151 variable (which in Fourier space reads TLf = ik/\k\f), 

^ p signature of the Landau damping. The specific form 

ft I . . b °f the expressions for the gyrotropic fourth rank cumu- 

F x = puu x + -l^sm a p l ± + cos a {p\ - P l ± )-^L hints and f r ±± , of the gyroviscous tensor II, of the 

non-gyrotropic contributions to the fourth-rank cumu- 

+ sina H xx + cos a IL^ J (A16) i a nt R r NG , of the transverse components of the fluxes of 

1 {3 / i i by parallel and perpendicular heat S x and S x as well as 

F y — puuy + — ^cosa (p|| — Px.) u p of the coefficients T^xp and function C\, are computed 

\ from the linear kinetic theory and given in They 

+ sina U xy + cos a U yz j (A17) i nvo i ve functions r and Ti where T u (b) is the product 

P , ■ ■ b °f cxp(— &) by the modified Bessel function I u (b), with 

F " = + 2 ( C ° S a * + C ° S a (P " - P1) W N kirl/2 = ^ sin 2 a. 



+ smaIL xz j (A18) 

-Fj = sin a cos a b x + ^ ^ sin a 

+ cosa (p|-pl)^) (A19) 

El = - cos a 6, + ^(cosa (pjj -pi)|fe) (A20) 

i*^ = cos a — cos a o 2 + — ^ cos a p ± 

+ cos«(pf-pl)|^). (A21) 
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FIG. 1: Contours of the e-pseudospectrum for a plasma 
with /3 = 10, k = 0.5, T±/Tj| = 1. Contours are plotted 
for logioe = —6.8, —6.6, . . . , —5.6. The dotted line is how 
the e-contour would appear if the operator were normal, for 
e = 0.00075. 
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FIG. 2: Evolution in time of the absolute value of the ampli- 
tude of the y component of the magnetic field, normalized to 
its initial value, for one particular choice of initial condition. 
The parameters used are: k = 1, T±/Tu — 1, /3 = 1 (top 
panel) and f3 = 5 (bottom panel). The curve in dashed line 
is the evolution e at , predicted by modal theory. 
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FIG. 3: Evolution in time of the absolute value of the ampli- 
tude of the y component of the magnetic field, normalized to 
its initial value, for one particular choice of initial condition. 
The parameters used are: T±/T« = 0.65, f) = 5, k = 1 (top 
panel), k = 5 (central panel), and k = 10 (bottom panel). 
The curve in dashed line is the evolution e at , predicted by 
modal theory. 




FIG. 4: (Color online) Average value of G(t) computed over 
10000 randomly generated initial perturbations, shown at four 
different times, in the (f3, l—T±/Tu) space. Logarithmic scale. 
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FIG. 5: (Color online) Average (left panel) and maximum 
(right panel) values of SB /Bo computed over 10000 randomly 
generated initial perturbations, shown at four different times, 
in the (/3, 1 — T±/T\\) space. Logarithmic scale 
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FIG. 6: Numerical abscissa r\ as a function of k for the fol- 
lowing parameters: (3 = 2 (black), 5 (red), 10 (green), and for 
Tx/T|| = 0.5 (dot-dashed), 1 (solid line), 2 (dashed). Note 
that the isotropic case is marked with solid line. 



